Inverse of Log CDF of Normal Distribution · Issue #13923 · scipy/scipy · GitHub

Is your feature request related to a problem? Please describe. We were interested in calculating z-scores associated to log p-values in the situation where log(p) can be very small. The naive solution of applying scipy.special.ndtri(numpy.exp(log_p)) fails when log_p < ~-745 due to underflow.

Describe the solution you'd like In the underlying C code for scipy.special.ndtri found here, when p is very small the approximation used for ndtri(p) is an expression in sqrt(-2*log(p)). This formula could then be applied directly to calculate the inverse of the Log CDF of the normal distribution for very small log(p). A function could be written which exponentiates log_p and then uses the existing ndtri implementation on the exponentiated value, except for when log_p is very small, in which case the rational approximation in sqrt(-2*log_p) can be applied directly. Preferably this could be done in C or Cython for performance reasons.

Describe alternatives you've considered

We've already implemented this approach in Python and applied it in our research. A code snippet can be found below. Would there be interest in adding such a solution to scipy through a C extension? If so, I'd be happy to submit a PR, though it may be several weeks before I have the free time to do so. In the mean time, if anyone else has an interest in adding this that would be great too.